---
title: "Lab 13 Code"
author: "Fabio Setti"
date: "Spring-2025"
output: html_document
---


## Today's Packages and Data 🤗

```{r}
install.packages("naniar")
install.packages("missMethods")
```

```{r}
library(naniar)
library(missMethods)
library(tidyverse)
theme_set(theme_classic(base_size = 16, 
                        base_family = 'serif'))
```

```{r}
dat <- rio::import("https://github.com/quinix45/PSYC-7804-Regression-Lab-Slides/raw/refs/heads/main/Slides%20Files/Data/employee.csv")
```


```{r}
#lav_summary() function

source("https://github.com/quinix45/PSYC-7804-Regression-Lab-Slides/raw/refs/heads/main/Slides%20Files/Additional%20files/lavaan_summary_lab12.R")
```

## Software Behavior With Missing Data?

```{r}

climate_miss <- c(rep(NA, 300), dat$climate[301:630])

# Quick way of counting missing values in a column
sum(is.na(climate_miss))
```

```{r}
reg_miss <- lm(empower ~ lmx + worksat + climate_miss, dat)
length(resid(reg_miss))
```

## Example of Bias: The Population

```{r}
set.seed(20345)

# generate some data that would show bias due to restriction of range

X <- rnorm(1000)
Y <- rnorm(1000, mean = ifelse(X > .7, .2*X, .4*X), sd = .4)

bias_dat <- data.frame("Graduate_GPA" = Y,
                       "GRE_Score" = X)
```

```{r}
cor(bias_dat$Graduate_GPA, bias_dat$GRE_Score)
```

```{r}
ggplot(bias_dat,aes(x = GRE_Score, y = Graduate_GPA))+
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
         labs(title = "Unbiased Sample of the Population") +
   theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16))
```

## Example of Bias: Biased Sample 


```{r}
cor(bias_dat[bias_dat$GRE_Score > 0,])[2]
```


```{r}
ggplot(bias_dat[bias_dat$GRE_Score > 0,] , aes(x = GRE_Score, y = Graduate_GPA))+
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
         labs(title = "Biased Sample of the Population",
            y= "GPA", 
            x = "Test Score") +
   theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16))
```


## Little's MCAR test


```{r}
# random missingness
library(missMethods)

set.seed(78984)
# 30% missing
dat_miss_MCAR <- delete_MCAR(dat, p = .3)
```


```{r}
# run test (need to run code above first)
mcar_test(dat_miss_MCAR)
```

```{r}
set.seed(467)
# Probabilities of missingness based on `worksat` 
mis_prob <- 1 - plogis(dat$worksat, location = mean(dat$worksat), scale = sd(dat$worksat))
dat_miss_MAR <- dat

# MAR process for other variables
for(i in 1:4){
  dat_miss_MAR[,colnames(dat)[-3][i]] <- ifelse(rbinom(nrow(dat), 1, mis_prob) == "1", NA, dat[,colnames(dat)[-3][i]])
}
```

```{r}
# run test (need to run code above first)
# `worksat` taken out because I left it complete on purpose
mcar_test(dat_miss_MAR[,-3])
```

## Non-Significant Little's Test?

```{r}

reg <- lm(worksat ~ empower + lmx + climate + cohesion, dat)

round(coef(reg), 3)

```


## MCAR VS MAR with `lm()`


```{r}
# MCAR sim


set.seed(78984)

MCAR_list <- c()

for (i in 1:2000) 
  {
  dat_miss_MCAR <- delete_MCAR(dat, p = .3, cols_mis = colnames(dat)[-3])  
  
  MCAR_list[[i]] <-  coef(lm(worksat ~ empower + lmx + climate + cohesion, dat_miss_MCAR))
  
  }

MCAR_res <- bind_rows(MCAR_list)

```

```{r}
round(colMeans(MCAR_res), 3)
```


```{r}

# MAR sim

set.seed(4627)


MAR_list <- list()
mis_prob <- 1 - plogis(dat$worksat, location = mean(dat$worksat), scale = sd(dat$worksat))  
  
for(i in 1:2000){

dat_miss_MAR <- dat

# MAR process for other variables
for(j in 1:4){
  dat_miss_MAR[,colnames(dat)[-3][j]] <- ifelse(rbinom(nrow(dat), 1, mis_prob) == "1", NA, dat[,colnames(dat)[-3][j]])
}

MAR_list[[i]] <-  coef(lm(worksat ~ empower + lmx + climate + cohesion, dat_miss_MAR))

}

MAR_res <- bind_rows(MAR_list)
```


```{r}
round(colMeans(MAR_res), 3)
```

## FIML with `lavaan`

```{r}
library(lavaan)
reg_mod <- "worksat ~ empower + lmx + climate + cohesion"

fit_miss <- sem(reg_mod, dat_miss_MAR, meanstructure = TRUE, 
                fixed.x = FALSE, missing = "fiml")

# parameterestimates(fit_miss)[16, 4] gets the intercept for this specific model (do not copy and paste this line and expect it to print the intercept for other models, it will not work; you need to look for it through the lavaan summary yourself usually)
round(c(parameterestimates(fit_miss)[16, 4], 
        lav_summary(fit_miss)[,4]), 3)
```



```{r}
# compare to original values

round(coef(reg), 3)
```

## Comparing `lm()` and `lavaan` With MAR data


```{r}
round(coef(reg), 3)
```


```{r}
set.seed(4627)


MAR_list_lm <- list()
mis_prob <- 1 - plogis(dat$worksat, location = mean(dat$worksat), scale = sd(dat$worksat))  
  
for(i in 1:2000){

dat_miss_MAR <- dat


# MAR process for other variables
for(j in 1:4){
  dat_miss_MAR[,colnames(dat)[-3][j]] <- ifelse(rbinom(nrow(dat), 1, mis_prob) == "1", NA, dat[,colnames(dat)[-3][j]])
}

MAR_list_lm[[i]] <-  coef(lm(worksat ~ empower + lmx + climate + cohesion, dat_miss_MAR))

}

MAR_res_lm <- bind_rows(MAR_list_lm)
```


```{r}
round(colMeans(MAR_res_lm),3)
```



```{r}
# this will take around 10 minutes to run (could be faster if I parallelized loops)
# make sure the lav_summary() function is loaded


library(lavaan)

set.seed(4627)

MAR_list_fiml <- list()
mis_prob <- 1 - plogis(dat$worksat, location = mean(dat$worksat), scale = sd(dat$worksat))  
lav_mod <- "worksat ~ empower + lmx + climate + cohesion"
  
for(i in 1:2000){

dat_miss_MAR <- dat

# MAR process for other variables
for(j in 1:4){
  dat_miss_MAR[,colnames(dat)[-3][j]] <- ifelse(rbinom(nrow(dat), 1, mis_prob) == "1", NA, dat[,colnames(dat)[-3][j]])
}


fit_miss <- lavaan::sem(lav_mod, dat_miss_MAR, meanstructure = TRUE, fixed.x = FALSE, missing = "fiml")

MAR_list_fiml[[i]] <-  c(parameterestimates(fit_miss)[16, 4], lav_summary(fit_miss)[,4])

print(paste("iter", i))

}

MAR_res_fiml <- matrix(unlist(MAR_list_fiml), ncol = 5, byrow = TRUE)
```

```{r}
round(colMeans(MAR_res_fiml), 3)
```



